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(~| , Abstract 

O 

1^ I We show that integrabihty of the BCS model extends beyond Richardson's model (where 

all Cooper pair scatterings have equal coupling) to that of the russian doll BCS model for 
which the couplings have a particular phase dependence that breaks time-reversal symme- 
Cu . try. This model is shown to be integrable using the quantum inverse scattering method, 

c/3 ' and the exact solution is obtained by means of the algebraic Bethe ansatz. The inverse 

problem of expressing local operators in terms of the global operators of the monodromy 

e matrix is solved. This result is used to find a determinant formulation of a correlation 

function for fluctuations in the Cooper pair occupation numbers. These results are used 
1-^ ' to undertake exact numerical analysis for small systems at half-filling. 

a' 

o 

O ; 1 Introduction 






Experimental work conducted by Ralph, Black and Tinkham [1,2] was successful in deter- 
^ I mining the discrete energy spectra of aluminium nanograins with a fixed number of electrons. 

Their results surprisingly showed remnants of superconducting behaviour. Significant dif- 
CN . ferences between the spectra of grains with even or odd number of electrons suggested that 

^^ I strong pairing correlations characteristic of superconductivity were still present, and quantum 

t:j- ' fluctuations significant. The Bardeen, Cooper and Schrieffer (BCS) model of superconductiv- 

^D ■ ity [3] can be applied in this situation, but the usual mean-field approach is not suitable for 

two reasons: firstly the number of electrons must be fixed, and secondly the strong quantum 
fluctuations destroy the validity of the assumption that operators may be averaged. 

Fortunately predictions from Richardson's model, for which all scatterings between pairs 

'^ i are assumed to have equal coupling, give an excellent match with the experimental results 

Q I discussed above [4]. In particular, this model is exactly solvable. The solution was originally 

O ' obtained and analysed in the 1960s by Richardson and Sherman in the context of nuclear 

physics [5,6]. The model has also been shown to be integrable [7], and by recasting this result in 

K^ \ the context of the quantum inverse scattering method, the solution was subsequently obtained 

^ • by means of the algebraic Bethe ansatz method [8,9]. The latter approach has the advantage 



cd 



C^ 



that it is possible to derive exact expressions for form factors and correlation functions [8, 
10, 11]. These considerations have already been used to produce several generalisations of 
Richardson's model [12-14], and to study other pairing Hamiltonians [15, 16]. Reviews on 
these topics can be found in [4, 11, 17]. Exact studies of these BCS models are also important 
for studying pairing correlations in systems of ultracold dilute Fermi gases [18]. 
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Here we study the russian doll BCS model, a one-parameter generalisation of Richardson's 
model. It differs from Richardson's model by the inclusion of phases (indexed by a parameter 
a) in the pair scattering couplings, which break time-reversal symmetry. Richardson's model 
is reproduced in the limit as the parameter a is taken to zero. The terminology 'russian 
doll' refers to the russian craft of making nested wooden dolls, known as matryoshka. Inside 
each doll one finds a smaller version of the doll. A russian doll renormalisation group flow 
is one which displays a cyclic nature rather than flowing to a fixed point. The russian doll 
BCS model was proposed in [19] as an example of a many-body system exhibiting a russian 
doll renormalisation group flow. In §2, the model is shown to be integrable via the quantum 
inverse scattering method and its exact spectrum is obtained via the algebraic Bethe ansatz 
method. Our derivation shows that the Hamiltonian is embedded in the transfer matrix as the 
second-order term in the expansion about an infinite spectral parameter. This construction 
is in contrast to the proof of integrability for Richardson's model given in [11] in that the 
quasi-classical limit is not taken. In this sense there is no connection between the russian doll 
BCS model and Gaudin magnets (cf. [17]). The effect of the parameter a on the spectrum of 
the model is discussed in §3. In particular we will show that a influences the degeneracies of 
the energy bands in the strong coupling limit. In §4 we solve the inverse problem. This result 
is used in §5 to obtain an exact formula for a superconducting pairing correlation function. 
Numerical analysis of our results show that the ground state of the model is qualitatively 
independent of a in terms of the fluctuations in Cooper pair number for the single particle 
energy levels. Increasing a does however lead to a quantitative suppression of the fluctuations 
across all levels. Our conclusions are given in §6. 

2 Integrability of the russian doll BCS model and its exact 
solution 

The physical properties of a metallic nanograin with pairing interactions are described by the 

BCS Hamiltonian [4] 

c c 

Hbcs = J2£jnj - Yl 9jA+A-(^j-Cj+- (1) 

Here, j=l, . . . , £ labels a shell of doubly degenerate single particle energy levels with energies 
£j, and Uj = Cjj^c.j^ + (yj_c,_ is the fermion number operator for level j. The operators 

c_|_, cj_|_ are the annihilation and creation operators for the fermions at level j. The labels ± 
refer to pairs of time-reversed states. 

An important aspect of the Hamiltonian (^ is the blocking effect. For any unpaired 
fermion at level j the action of the pairing interaction is zero since only paired fermions are 
scattered. This means that the Hilbert space can be decoupled into a product of paired and 
unpaired fermion states in which the action of the Hamiltonian on the space for the unpaired 
fermions is automatically diagonal in the natural basis. In view of the blocking effect, it is 
convenient to introduce hard-core boson operators h- = c-c-, , Ir- = (yj,c}j_ which satisfy the 

{h]f = 0, [5^., bl] = 8^,{l - 2b]b^), [b^, 6J = [h], &t] = 

on the space excluding single particle states. We also set N- = b-b, to be the Cooper pair 
number operator for the j'th level. In this setting the hard-core boson operators realise the 



(2) 



su{2) algebra in the pseudo-spin representation, through the identification 

S- = b, 5+ = 6^ S' = ^i2N-I). 

The pseudo-spin operators satisfy the following commutation relations: 

[5^ S^] = ±S^, [S+, S-] = 2S'. (3) 

It is now well known that the Hamiltonian ^ can be solved exactly in the case when 
9jk = 5' is constant for all I < j,k < C. This is referred to as Richardson's model. The 
solution was obtained in [5]. It was subsequently shown that the model is also integrable 
through an explicit construction for a set of conserved operators [7]. Both of these results 
can be obtained from the formulation of the model through the quantum inverse scattering 
method [8,9]. Here we will adopt this approach to establish that there exists a more general 
manifold of integrability. We will show that the following Hamiltonian 
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2^e,Ar,-5 5](e-6i6^.+e— 6tft^ 



(4) 



i=i 



j<k 



is integrable for arbitrary values of a. This model includes Richardson's model in the limit 
a — »■ 0, up to the addition of a multiple of the u{l) operator N = X],=i Nj for the total number 
of Cooper pairs. It coincides with the russian doll BCS model recently studied in [19]. 

To establish the integrability of (jlj), we begin with the i?-matrix R(u) that satisfies the 
Yang-Baxter equation 



i?i2(n - v)Ri3{u)R23iv) 
and has the form [20,21] 



R23{v)Ri3{u)Ri2{u 



(5) 



R{u) 



1 



-{ul ®I + iriP) 
u + irj 

/ 1 \ 

h{u) c{u) 

c{u) h{u) 

\ 1/ 



(6) 



where h[u) = u/{u + ir]), c{u) = 7]/{u + irj) and 7] is taken to be an arbitrary real parameter. 
Above, P is the permutation operator. It satisfies 

P{x ® y) = y (8) X 

for all vectors x and y. An important property of the i?-matrix that will be exploited later is 
that it satisfies the unitarity condition 



R{u)R{-u) = 
The i?-matrix is (7/(2)-invariant in that 

[R{u),g^Q] 



1(^1. 



(7) 



(8) 



where q is any 2x2 matrix. We introduce the Lax operator L(u) in terms of the su{2) Lie 
algebra with generators S^ and S^ [22,23], 

^^ u + irjy ir]S+ ul + ir]{I/2 - S') J ' ^' 

subject to the commutation relations (jSJ- Then the following holds: 

Ri2{u - v)Li{u)L2{v) = L2{v)Li{u)Ruiu - v). (10) 

When the su{2) algebra takes the spin 1/2 representation Q the resulting L-operator is 
equivalent to that given by the i?-matrix ©. 

We take g = exp(— iao") with a = diag(l, —1) and use © to construct the monodromy 
matrix 

T{u) = QoLociu- ec)---Lo2iu- e2)Loi{u-ei) 

' A{u) B{u) 
C{u) D{u) 

Above, labels the auxiliary space, while the tensor components of the physical space are 
labelled 1, . . . , i2. Thus the entries A{u), B{u), C{u), D{u) of T{u) are elements of the £-fold 
tensor algebra of su{2). As a consequence of (jH| and (|lfl|) we have the relation 

Ri2{u - v)Ti{u)T2{v) = T2{v)Ti{u)Ri2{u - v). (11) 

We next define the transfer matrix 

t{u) = troT(ti) 

= A{u) + D{u). 

Here, trg denotes the trace taken over the auxiliary space. It can be shown using Hll|) that 
the transfer matrix t{u) generates a one-parameter family of commuting operators, viz. 

[t{u),t{v)]=Q \/u,veC. 

Next we specialise the representation of the su{2) algebra to that given by 0. Following the 
method of the algebraic Bethe ansatz (see for example [11]) we can deduce that the transfer 
matrix eigenvalues are 

» / N / . N TT u- ^k j-r u-Vj + 3iri/2 , . ^ tt "" - ^i - ^^/2 

A{u) = exp(-za) — ^ ^ + exp(za) ^ V 

f^J- u — £k + if] -Hr u — Vj + irj/2 -H- u — Vj + ir]/2 

subject to the constraints of the Bethe ansatz equations 

exp^ la) jj ^. _ £^ + i^/2 ii y._y.^iri- ^ ' 



It is convenient to define an equivalent transfer matrix 

n — Cfc + iry 



t(n) = n 



L u -Efc + ir]/2 ^ 



tiu) 



whicli has eigenvalues 



Hu) = ~a{u) \\ ^--^+^^^1^ + d{u) n " - "^- - ^^/' 



where 



c 



U-Sk 



~a{u) = expHa)n^_,, + ,^/2 

-c , ■ 

u- Ek+tr] 



d{u) = exp(ia) IT 



L u-efc + i'n/2' 



(13) 
(14) 



Analogously we define 



i(n) = n 



u — Ek + if] 



A{u) 



^^J-^u-ek + iv/2^ 

and similarly for B{u), C{u) and D{u). 

A series expansion for t(n) yields a set of mutually commuting operators. Here we choose 
to expand about the point at infinity: 



i{u) = Y.~^ 



{k)^-k 



fc=0 



such that 



It is straightforward to deduce the leading terms: 

(■ 2 —2 
/ + riu~^ tan(a)(7V - C/2) - !!L_^ tan(a)(iV - C/2) 

c 



+ 7JU tan(a 
1 



E^^^-^E^ 



i=l 



2 -2 

77 n 



iN-C/2Y 



i=l 

C 1 

+ 



8 2 cos (a) 



X;(e^"6l6,+e— 6]6, 



Expanding A(ii) in powers of n ^ we find 



j<k 



,2„,-2 



(15) 



A(u) ~ 2cos(a) l + r/u"Han(a)(iV-£/2) tan(a)(iV - /:/2) 
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Comparing this expression with the expansion (|15j) we immediately see that the exact solution 
for the energies of the Hamiltonian (0J) with g = r//sin(a) is 

N 

E = 2'^Vj +gN cos(a), 

where the {vj } are solutions to p2p . Making the change of variable 

a ^ r]a 
and taking the limit r/ — j- this reproduces the known exact solution [5], given that 

As a set of conserved operators for the Hamiltonian @ we take {t{£k) '■ k = 1, . . . , C}, which 
naturally generalise those of [7] , and arise in the solution of the inverse problem below. These 
operators clearly satisfy 

[Hrd, t{ek)] = [t{ei), t{ek)] =0, yk,l = l,...,£. 

Hence the above construction shows that the russian doll BCS model is both integrable and 
exactly solvable. 

We note that the unitary time-reversal transformation of @ is equivalent to the change 
of variable a -^ —a with g and each Ej fixed. It is easily seen that the change a — > —a and 
ry — > — 7/ leaves g, E and each of the Bethe ansatz equations (|12jl invariant. We mention also 
that in the above the single particle energies £j are entirely arbitrary. In particular, it is not 
necessary to have Ej < Ek for j < k. 

3 Energy spectrum 

To study the behaviour of the spectrum of Q as a varies between and tt/2, we solved the 
Bethe ansatz equations (J12j) for the picket fence model in which the energy levels are equally 
spaced, that is Ej = j — (£ + 1)/2. With this convention the Fermi level is zero for a system at 
half- filling. The blocking effect simplifies the calculation of the energy of the excited states, 
since a state consisting of say N' Cooper pairs and two free electrons blocking levels m and 
n has energy 

N' 
E = 2'^Vj + gN' cos(q) + Em + Sn 

J=l 

where the {v} = {vi, . . . ,V]\fi} satisfy the following Bethe ansatz equations 
exp(-2za) IT ^^i^It^Ml = fl !!l^^l^^ , 

For a fixed number of particles the excitations of the model Q are classified according 
to the initial distribution of the Cooper pairs and blocked energy levels at g = 0, thereby 



suggesting the initial configuration of the Bethe ansatz roots. For example, a solution of the 
Bethe ansatz equations for the ground state at half-filling is found using the initial conditions 
Vj = Ej , j = 1, . . . , C/2, whereas an excited state formed by promoting a Cooper pair above 
the Fermi level is obtained from the initial conditions Vj = £j , j = 1, • • • , {C—2)/2 and 

Vc/2 = e(£+2)/2- 

We solved the Bethe ansatz equations iteratively starting with each possible configuration 
of the roots at 5 = and a = 0. Using this method we calculated the excited state spectrum 
for a model with three Cooper pairs at half- filling (£ = 6), a total of 141 states if the particle- 
hole symmetry is not taken into account. Figure ^ shows the energy levels as a function of 
the coupling g for three different values of a. For small a and large coupling g, the energy 
levels form narrow, well-separated bands in agreement with the results of Richardson's model 
(e.g. see [24]). As a increases the energy levels move, the original bands split apart, and new 
bands are formed. Throughout, the gap between the ground state energy and the first group 
of excited states remains, though it decreases as a increases. There is also a set of states for 
which the energy of each state does not depend on a. These are the states that consist of 
one Cooper pair and all levels apart from levels I and m are blocked. The single Bethe ansatz 
equation has two solutions: 

2v± = ei+em± [g^ + {ei - emf] - gcos{a). 



Consequently the energy for this state is 

E±=ei + em± [g^ + (ei - £mfj + ^ £i, 

where B indexes the set of blocked energy levels. 

Despite the fact that the Bethe ansatz equations have to be solved numerically to obtain 
the root positions for nonzero values of g, Gaudin [25] and Sierra et al. [26] have made some 
conjectures as to the behaviour of the BCS roots at very large g for Richardson's model. 
Numerical studies indicate that as g increases some of the roots form complex conjugate pairs 
in a complicated pattern which may involve roots forming pairs, splitting apart, and forming 
new pairs, or complex conjugate pairs separating and becoming independent real roots again. 
The root behaviour of the ground state of a model with an even number of Cooper pairs is 
very simple: all roots form complex conjugate pairs and become infinite at large coupling. 
For an odd number of Cooper pairs the roots behave in the same way apart from the root 
initially closest to the Fermi level, which remains finite. Generally most roots tend to infinity 
at very large g whether or not they have formed a complex conjugate pair. The number of 
the remaining finite roots can be thought of as the number of elementary excitations of the 
excited state. There is a formula due to Gaudin [25] that predicts the number of states with 
a particular number of finite roots at 5^ = 00, and an algorithm due to Sierra et al. [26] that 
matches the initial configuration of roots (or distribution of Cooper pairs over the available 
energy levels) to a final state. 

It is clear from the splitting of the energy bands shown in Figure ^ that Gaudin's formula 
and the algorithm due to Sierra et al. are not valid for nonzero a. In fact the russian doll 
BCS model exhibits the same general root behaviour as Richardson's model (as shown in 
Figures IS^-c) but many more roots tend to infinity as g ^ 00. In each of the figures the large 
g configuration of Bethe ansatz roots consists of two real roots and one complex conjugate 
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Figure 1: The behaviour of the energy levels as a is increased. The results shown are for 
the picket fence model of £ = 6 single particle energy levels at half-filling. As the coupling g 
increases, the energy levels coalesce into bands. Variations in a lead to different degeneracies 
of the bands. Note that in each case the ground state is non-degenerate and there is no level 
crossing between the ground and first excited energy level as g increases. 



pair of roots. For a = 0.01, Figure|2^ shows there are two finite real roots, in agreement with 
the algorithm which predicts two elementary excitations [26]. At a close to 7r/2 (see Figure 
Et), it appears that a single real root remains finite, but upon close examination at g' ~ 100, 
the root is seen to be very slowly diverging. Thus there are no longer any finite roots in this 
particular case. 

4 Solution of the inverse problem 

In preparation for the next section we would like to be able to express any local operator in 
terms of the elements of the monodromy matrix. This is known as the inverse problem, and 
its solution will enable the form factors of the local operators to be determined. The results of 
this section are the generalisation of those in [27] to the case where the operator g is included 
in the definition of the monodromy and transfer matrices (cf. [8]). 

Hereafter we set Rjk = Rjk{£j — £k) for notational ease. First we note that 



T{ej) - Rj{j-i) ■ ■ ■ RjiPojQjRjC ■ ■ ■ Rj(j+i) 
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Figure 2: The evolution of the real and imaginary parts of the roots {v} with increasing g. 
The three plots shown are for fixed values of a. In each case the roots correspond to the excited 
state of the C = 8 model with A^ = 4 Cooper pairs and initial conditions {v} = {ei, £3, £5, eg}- 



or equivalently 

Poj = Rji ■ ■ ■ ^7(i-i)^(^i)^7{j+i) ■ ■ ■ ^JcQj ■ (16) 

It is readily deduced that 

t{ej) = Rj(j-i) ■ ■ ■ RjiQjRjc ■ ■ ■ Rj{j+i)- (17) 

At this point we appeal to the unitarity condition which allows us to write 

t{ej) = Rlj_i)j ■ ■ ■ Rij QjRjc ■ ■ ■ Rjij+i)- 
Defining 

Tj-fc = RlkRlk • • • Rjk, 

we then have 

which leads to the formula 

k f ^ \ 

n *(^i) = n Si ^kcTk{c-i) ■ ■ ■ rfc(fc+i). (18) 

j=i \i=i / 

Equation ((TH|) is proved by induction, noting first that 

^(^l) = 01^1£-Rl(£-l) • • --^12 

= 0iri£ri(£_i) • ••ri2- 



Employing equation (|18j) we find 
Combining (|16|) and (|17|) we have 

-Pofc = Rkl ' ' ' ^A:(A;-l)^(^'«)* {^k)Rk{k-l) ' ' ' Rkl 

from wliich we obtain tlie following expression 

T{ek)t~ (efc) = '^ (k_i^k^Qk^ {k-^)k- 
Thus using 1)19^ we have 

Pok= \\\t{e,)Y{eJ\{t-\e,)\. (20) 

Finally, since 
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formula ()2U() leads to expressions for the local operators 6^, b\, Nj^ terms of the global operators 
A{u), B{u), C{u) and D{u) of the algebraic Bethe ansatz. 

5 Pair occupation fluctuations 

In this section we use the notation (x) = ({f}|x|{^}) / ({^}|{i'}) for any operator x where 
\{v}) is a Bethe state associated with the roots {v} = {vi, . . . ,viy}- For the BCS model it is 
of interest to study the behaviour of the correlation function C^ describing the fluctuations 
in the Cooper pair occupation numbers for each level k: 

Cl = {Nl)-{N,f 
= {I-Nk){Nk). 

For Richardson's model (a = 0), this correlation function has been studied in [4, 10]. 

Solving the inverse scattering problem allows the correlator to be expressed in terms of 
the form factor of the global operator D{u) via 

{{v}\l-N,\{v}) = {{v}\D{e,)t-He,)\{v}). 

Fortunately the Slavnov formula [28] (see also [11,22,29]) for the scalar products of Bethe 
states leads to explicit determinant representation for the form factors of the operators 
A{u), B(u), C{u) and D{u). We follow [11], in which the Slavnov formula is described for 
t(u — ir]/2) together with form factors for the elements of the associated monodromy matrix 
(up to a rescaling of rf) . We define it here in terms of the shifted functions 

a{u) = a(u — irj/2)^ d(u) = d(u — irj/2) 

10 



where d{u) and d{u) are given by (|T3|) and (fTH) respectively. 

Below both sets of parameters {wi} and {vi} are assumed to be solutions of the Bethe 
ansatz equations (jT^ . The scalar products of the Bethe states are obtained from 

({^}iw) = Ywi hw~t ;' 

where i*" is an A^ x iV matrix with entries 

N N 



^ indiwi) 



{Vj - Wi) 



JJ(uj -Wk+ irf) - d{vj) Wivj -Wk- irf) 



unless {w} = {v}, in which case the diagonal entries become 

a{vi) 



d{v.) Y. 



Vi-Vl+IT] 
l^i k=l 

/ N 






Vi-Vl-IT] 
l^i k=l 

+ d{vi) a'{vi) W_{vi -Vk + irf) + d'{vi) W_{vi -Vk-irf)]. 

\ k=l k=l J 



(21) 



where the prime denotes the derivative. Note that (|21|) corrects typographical errors appearing 
in equation (46) of [11]. The matrix elements of the operator D{u) were derived in [11] and 
read 

/ N 



{{w}\D{u)\{v}) = —j^- 



n :^"!^:-!(' idet(F+Q(.)) 



where 
and 

Qij{u) 



N 



= TT u-Wk + ^r^/2 
ill u-Vk + ivl2 



irid{wi)d{vj) 



{u — Wi + irj/2){u — Wi — iri/2) 



N 



d{u) -pr u — Wk + 3iri/2 
d{u) y. u-Wk- irj/2 



N 



Ylivj -vi-ii]). 



In the limit u ^ £k the term within the square brackets becomes equal to unity since a(efc) = 0, 
and the expression for the form factor of the local operator Nk between identical Bethe 
eigenstates is also very simple: 



(Nk) 



^ detjF + QJEk)) 



detF 



Figure 131 shows the behaviour of Ck for the ground state of the C = 6 system at half-filling, 
for three values of the coupling strength g as the parameter a is increased from just above 
zero to just below 7r/2. 
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Figure 3: The correlation function C^ for various choices of a and the couphng strength g. 
The results shown are for the ground state of the i2 = 6 model at half-filling. In each case the 
distribution of Cfc peaks about the Fermi level. For fixed g, increasing values of a suppress 
the magnitude of C^ across all single particle energies e^. 



6 Conclusion 

We have demonstrated the integrability of the russian doll BCS model [19], which generalises 
Richardson's model by the inclusion of phases in the pair scattering couplings, parameterised 
by a variable a. We have applied the algebraic Bethe ansatz to find the exact solution of 
the model, and have also solved the inverse problem for the computation of form factors and 
correlation functions. In particular, an explicit expression for the one-point correlator char- 
acterising Cooper pair occupation fluctuations was obtained. Analysis of this result indicates 
that the ground state structure is qualitatively independent of q, and increasing a suppresses 
the fluctuations across all single particle energy levels. On the other hand, our analysis of 
the energy spectrum shows that the degeneracies of the energy bands in the strong coupling 
limit is dependent on a. An open problem is to determine the extent to which a affects the 
thermodynamic properties of the model. The quantum transfer matrix method and the solu- 
tion of the associated nonlinear integral equations constitute powerful techniques for exactly 
calculating the free energy of many integrable systems. These techniques have successfully 
been applied to the cases of the XXZ model [30] and an integrable spin ladder [31]. As we 
have shown that the Hamiltonian is embedded in the transfer matrix as the second-order term 
in the series expansion about infinite spectral parameter, we hope that our results are a first 
step in formulating the quantum transfer matrix method for the russian doll BCS model. 
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